clear all

cd "$revdta"
set more off


***********PARTICIPATION MOMENTS BEFORE AND AFTER AGE 45 ****************************************************************************************************************
scalar def lowlim		=0.85
scalar def highlim		=5

local type $type
local typename $typename

use ".\Full_Panel_2b.dta", clear
keep if sample_analytic==1 & educ <= 2
bysort newid: egen firstage = min(age)
bysort newid: egen firsthealth = max(DS * (age==firstage))

*local var firsthealth
*local var chd_maxhealth
local var

*I lose most people before 1996 from child health, another 30% from restricting 
*to in-sample heads
*total: 2k households

g sp_age = agew
g sp_age2 = agew^2/100


gen sp_dis = sp_DS==1 | sp_DS==2

gen part_full = lw!=. & hours >=1500 
gen part_part = lw !=. & hours < 1500 & hours >=500

gen sp_part_full = sp_lw !=. & sp_hours >=1500
gen sp_part_part = sp_lw != . & sp_hours < 1500 & sp_hours >=500


replace sp_DS = sp_DS >=1
gen     age_band=1 if age<=45
replace age_band=2 if age>45

gen     sp_age_band=1 if agew<=45
replace sp_age_band=2 if agew>45 & agew != .
tab age_band,gen(ad)


merge m:1 newid using "./intermediary/sp_temp_types"
keep if _merge ==3
ren type_fix type_lp
drop lw_res


ren $type type_fix

tab type_fix, gen(prodtype_)
g ageh = age*(mod==0)*(sev==0)
g sp_ageh = sp_age*(sp_dis==0)*(sev==0)
g agemod = age*mod
g sp_agemod = sp_age*sp_dis
g agesev= age*sev

g age2h = age2*(mod==0)*(sev==0)
g sp_age2h = sp_age2*(sp_dis==0)*(sev==0)
g age2mod = age2*mod
g sp_age2mod = sp_age2*sp_dis
g age2sev= age2*sev

g old = age >= 45

g agehold = age 
replace agehold = 0 if age >= 45 | sev == 1 | mod == 1
g age2hold = age2 
replace age2hold = 0 if age >= 45 | sev == 1 | mod == 1

g sp_agehold = sp_age
replace sp_agehold = 0 if sp_age >= 45 | sp_dis == 1
g sp_age2hold = sp_age2 
replace sp_age2hold = 0 if sp_age >= 45 | sp_dis == 1

g agemodold = age 
replace agemodold = 0 if age >= 45 | mod == 0
g age2modold = age2 
replace age2modold = 0 if age >= 45 | mod == 0 

g agesevold = age 
replace agesevold = 0 if age >= 45 | sev == 0
g age2sevold = age2 
replace age2sevold = 0 if age >= 45 | sev == 0 

g sp_agemodold = sp_age 
replace sp_agemodold = 0 if sp_age >= 45 | sp_dis == 0
g sp_age2modold = sp_age2 
replace sp_age2modold = 0 if sp_age >= 45 | sp_dis == 0 


g hold = age >= 45 & sev == 0 & mod == 0
g sev_old = 0
replace sev_old = sev if age >= 45
g mod_old = 0
replace mod_old = mod if age >= 45

g sp_old = sp_age >= 45

gen sp_sev = sp_DS==2
gen sp_mod = sp_DS==1
g sp_hold = sp_age >= 45 & sp_dis==0
g sp_mod_old = 0
replace sp_mod_old = sp_dis if age >= 45


g vold = age >= 50
g sev_vold = 0
replace sev_vold = sev if age >= 50
g mod_vold = 0
replace mod_vold = mod if age >= 50

g vvold = age >= 55
g diffsev_vvold = 0
replace diffsev_vvold = diffsev if age >= 55
g diffmod_vvold = 0
replace diffmod_vvold = diffmod if age >= 55
g sev_vvold = 0
replace sev_vvold = sev if age >= 55
g mod_vvold = 0
replace mod_vvold = mod if age >= 55

drop yrd* age2
tab year,g(yrd)
g age2=age^2/100	
***** Wage regressions ****
reg lw prodtype_* sev mod age age2 old sev_old mod_old married yrd* , nocons
predict lw_res, residuals

eststo hwreg

/*SPOUSE - Drop those with outlier wage growth, below 1/2 of mw, or with less than 3 obs, or outside relevant age range*/

sort sp_newid year
qby sp_newid:gen grly=(exp(sp_lw)/exp(sp_lw[_n-1]))-1
replace grly=. if grly==-1
gen suspect=0
qui su grly,d

replace suspect=grly<-lowlim|(grly>highlim & grly!=.)	
egen ss=sum(suspect),by(newid)				/*# of outlier wage growth records*/
egen sw=sum(exp(sp_lw)<0.5*minwage),by(newid)	/*#of wages below state-level minimum wage*/
egen numy=sum(newid!=.),by(newid)			/*# of obs*/

gen sp_todrop=sw>0|numy<3|ss>0				/*Drop those with outlier wage growth, below 1/2 of mw, or with less than 3 obs*/
sort newid year

replace sp_lw = . if sp_todrop==1
replace sp_lw = . if agew < 23 | agew > 62

reg sp_lw prodtype_*  sp_mod sp_age sp_age2 sp_old sp_mod_old yrd1-yrd21, nocons

predict sp_lw_res, residuals

eststo swreg

***** Wage variance ***** 
qby newid: g laglw_res=lw_res[_n-1]
qby newid: g lagsp_lw_res=lw_res[_n-1]

gen lw_res2 = lw_res*lw_res
gen sp_lw_res2 = sp_lw_res*sp_lw_res
gen lwcov = lw_res*sp_lw_res

foreach x in "" "sp_" {
gen `x'lwlaglw_res = `x'lw_res*lag`x'lw_res
sum `x'lw_res if lag`x'lw_res != .
local mean1 = `r(mean)'
sum lag`x'lw_res if `x'lw_res != .
local mean2 = `r(mean)'
replace `x'lwlaglw_res = `x'lwlaglw_res - `mean1'*`mean2'
}

reg lw_res2
eststo wagevar

reg lwlaglw_res
eststo wagecov

reg sp_lw_res2
eststo spwagevar

reg sp_lwlaglw_res
eststo spwagecov

reg lwcov
eststo wagespcov



***** Average participation by age ****
matrix b=0
matrix V=0

foreach suffix in full {
         local i=0
	   while `i' < 3 {
		   local j=1	
		   while `j' < 3 {
				reg part_`suffix' if DS==`i' & age_band==`j' 
				eststo `suffix'_`i'_`j'
	         local j=`j'+1
                 }
         local i=`i'+1
         }
}
				 

foreach suffix in full part {
	forvalues m=1(1)1{
         local i=0
	   while `i' < 2 {
		   local j=1	
		   while `j' < 3 {
				reg sp_part_`suffix' if sp_DS==`i' & sp_age_band==`j' & married==`m'
				eststo sp_`suffix'_`i'_`j'
				local j=`j'+1
                 }
         local i=`i'+1
                 }
	}
}
			


********Stocks on DI by health, marital status, and age
sort newid age
bysort newid: gen worstDS = DS[1]
bysort newid: replace worstDS = max(DS[_n],worstDS[_n-1]) if _n > 1

		 local m = 0
		 while `m' < 2 {
         local i=0
	   while `i' < 2 {
		   local j=0	
		   while `j' < 3 {
				reg DI if DS==`j' & old==`i' & married==`m'
				eststo stockDI_m`m'_`i'_`j'
	         local j=`j'+1
                 }
         local i=`i'+1
                 }
				 local m = `m' + 1
				 }


**********COMPOSITION OF DI BENEFICIARIES BY HEALTH
tab DS, gen(DS)
		   local j=3	
		   while `j' > 0 {
		    local i=0
			while `i' < 2 {
				reg DS`j' if DI==1 & old==`i' 
				eststo compDI_`i'_`j'
			local i=`i'+1
                 }
	       local j=`j'-1
                 }


*********************************************







***************FLOWS ONTO DI BY MARITAL STATUS
xtset newid year
g lag1DI=l1.DI
g lag2DI=l2.DI
g lag1DS=l1.DS
g lag2DS=l2.DS

         local i=0
	   while `i' < 2 {
		   local j=0	
		   while `j' < 3 {
				reg DI if DS==`j' & ((lag2DI==0 & lag1DI == .) | lag1DI == 0) & old==`i'
				eststo flowDI_`i'_`j'
	         local j=`j'+1
                 }
         local i=`i'+1
                 }
				 
				 
				 

gen lc=ln(spending)

replace lc = . if year<=1996
save temp,replace

gen 	scale=1+0.5*(fsize-kids-1)+0.3*kids		/*OECD modified scale*/

					
replace lc=lc-ln(scale)
replace lc_impute=lc_impute-ln(scale)

gen mod_DI=mod*DI
gen sev_DI=sev*DI

gen mod_work = mod*part_full
gen sev_work = sev*part_full


gen DSDI0 =DS==0 &DI==1
gen DSDI1 =DS==1 &DI==1
gen DSDI2 = DS==2 &DI==1

 gen part_full_married = part_full*married
 gen part_full_single= part_full*(1-married)
 gen part_part_married = part_part*married
 gen part_part_single= part_part*(1-married)
 gen mod_single = mod*(1-married)
 gen sev_single = sev*(1-married)
 gen mod_married = mod*(married)
 gen sev_married = sev*(married)
 gen DSDI2_single = DSDI2 * (1-married)
 gen DSDI1_single = DSDI1 * (1-married)
 gen DSDI0_single = DSDI0 * (1-married)
 gen DSDI2_married = DSDI2 * (married)
 gen DSDI1_married = DSDI1 * (married)
 gen DSDI0_married = DSDI0 * (married)
 gen age_single = age*(1-married)
 gen age2_single = age2*(1-married)
 gen age_married = age*(married)
 gen age2_married = age2*(married)
 
 gen part_full_mod = part_full*mod
 gen part_full_sev= part_full*sev
 gen part_part_mod = part_part*mod
 gen part_part_sev= part_part*sev
 
 gen sp_part_full_dis = sp_part_full*sp_dis
 gen sp_part_part_dis= sp_part_part*sp_dis
 
  gen part_full_DS= part_full*DS
 gen part_part_DS= part_part*DS

  reg lc DS part_full part_part DI married part_full_DS part_part_DS sp_dis sp_part_full sp_part_part sp_part_full_dis sp_part_part_dis age age2 i.year
 eststo cons
 
 ****** Event study of health shock *********
gen baseyear = .
replace baseyear = -1 if year < 1996
replace baseyear = -2 if year >= 1996
 sort newid year
 gen disgroup =.
replace disgroup = 0 if sev==0 & mod == 0
replace disgroup = 1 if sev== 0 & mod == 1
replace disgroup = 2 if sev==1


bysort newid: egen onset_age_hsev = min(age+ 99999*(disgroup!=2)+99999*(sp_DS!=0))
replace onset_age_hsev = . if onset_age_hsev >= 9999
bysort newid: egen onset_age = min(age+ 99999*(disgroup==0))
replace onset_age  = . if onset_age >= 9999
bysort newid: egen anyonset_age = min(age+ 99999*((disgroup)==0))
replace anyonset_age = . if anyonset_age >= 9999

replace onset_age = . if onset_age>anyonset_age | onset_age == .

*Issue with these guys below is I don't know age of onset.
bysort newid: egen time1=max(age==onset_age-1)
bysort newid: egen time2=max(age==onset_age-2)
bysort newid: egen onset_year = max(year*(onset_age==age))
replace onset_age = . if time1==0 & onset_year <=1996
replace onset_age = . if time2==0 & onset_year >1996

sort newid year
gen emp = hours >=1500 & hours!=. 
bysort newid (year): gen empyet = sum(emp )
replace empyet = empyet>=1
bysort newid: egen empyet2=max(empyet== 1 & (age == onset_age + baseyear))
drop empyet
ren empyet2 empyet

 reg anyDI if age - onset_age >= 0  & age - onset_age <= 4  & onset_age_hsev < onset_age+5 & empyet==1
 eststo eventsev
 
  reg anyDI if age >= onset_age & age <= onset_age + 4   & empyet==1
  eststo eventmod

 suest hwreg swreg wagevar wagecov spwagevar spwagecov wagespcov cons full_* sp_full_* compDI_?_3 compDI_?_2 compDI_?_1 flowDI* eventmod eventsev,vce(cluster newid)
  
 mat def V = e(V)

local newnames
 local names: roweq V
local rownames: rownames V
 local rows = rowsof(V)
 local pos = 0
 tokenize `rownames'
 foreach name1 of local names {
		if(regexm("`name1'","lnvar") ==0 & ("`1'"!="_cons" | "`name1'" != "cons_mean")  &  regexm("`1'","part_part" )==0 & regexm("`1'","year" )==0 & regexm("`1'","yrd" )==0)    {
		local pos = `pos'+1 
		 local newnames `newnames' `name1'_`1'
	}
	macro shift
 }
 


  matrix finalV = J(`pos',`pos',0)
  local i = 0
  local pos1=1
 foreach name1 of local names {
 local varname1: word `pos1' of `rownames'
 if(regexm("`name1'","lnvar") ==0 & ("`varname1'"!="_cons" | "`name1'" != "cons_mean")  &  regexm("`varname1'","part_part" )==0 & regexm("`varname1'","year" )==0 & regexm("`varname1'","yrd" )==0) 	{
 local i = `i'+1
 local j = 1
 local pos2=1
	foreach name2 of local names {
	 local varname2: word `pos2' of `rownames'
		if(regexm("`name2'","lnvar") ==0 & ("`varname2'"!="_cons" | "`name2'" != "cons_mean")  &  regexm("`varname2'","part_part" )==0 & regexm("`varname2'","year" )==0 & regexm("`varname2'","yrd" )==0) 	{
		matrix finalV[`i',`j'] = V[`pos1',`pos2']
		local j = `j'+1
		}
		local pos2=`pos2'+1
	}
}
	  local pos1 = `pos1'+1
}
matrix rownames finalV = `newnames'
matrix colnames finalV = `newnames'
mat2txt, matrix(finalV) saving("$out/varcovmat_simplified") replace
